rm(list=ls(all=TRUE))
library(zoo)
library(plm)
library(reshape2)
library(data.table)
library(xlsx)
library(xtable)
library(DataCombine)
library(plyr)
library(vioplot)
library(ggplot2)
library(scales)
library(grid)
library(memisc)
library(sem)
setwd("~/Dropbox/LagIdentification/jop_final_version/replication")

#run simulation to replicate Figures 9 and A-4
simulation<-function(phi,rho,alpha,beta,gamma) {
    # panel properties
    #phi <- 1.5
    #rho <- 1.5
    #alpha <- 5
    #beta <- 1
    #gamma <- 1
    units<- 100
    time<- 50
    n<-units*time
    id <- seq(1:units)
    
    # generate fixed effects
    fe <- rnorm(units,0,1)
    
    # merge fes and create panel
    frame <- data.frame(id, fe)
    panel<-expand.grid(seq(1:time),seq(1:units))
    names(panel)<-c("t","id")
    total <- merge(frame,panel,by=c("id"), all=T)
    
    # function to generate individual panels of u = phi*u_t-1 + nu, nu~N(0,1)
    # and x = rho*x_t-1 + alpha*u + eta, eta~N(0,1)
    gen.panel<-function(nothing){
        z<-rnorm(time) #this is an instrument
        e2<-rnorm(time,0,1) #e1<-rnorm(time,0,1) # white noise
        e1<-vector(length=time) #e2<-vector(length=time)
        nu<-rnorm(time) # white noise
        e1[1]<-nu[1] #e2[1]<-nu[1] # random first value
        for(i in 2:length(nu)) {
          e1[i]  <- phi*e1[i-1] + nu[i]
        }
        #for(i in 2:length(nu)) {
        #    e2[i]  <- phi*e2[i-1] + nu[i]
        #}
        x<-vector(length=time)
        x[1]<-e2[1]
        for(i in 2:length(e2)) {
            x[i]  <- (rho*x[i-1] + alpha*e1[i] + gamma*z[i] + e2[i])/(1 - alpha*beta)
        }
        lagx<-vector(length=time)
        lagx[1]<- NA
        for(i in 2:length(e2)) {
            lagx[i]  <- x[i-1]
        }
        y<-vector(length=time)
        y[1]<-e1[1]
        for(i in 1:length(e1)) {
            y[i]  <- beta*x[i] + e1[i]
            #y[i]  <- (rho*x[i-1] +  + beta*gamma*z[i] + beta*e2[i] + e1[i])/(1 - alpha*beta)
        }
        lagy<-vector(length=time)
        lagy[1]<- NA # random first value
        for(i in 2:length(e1)) {
            lagy[i]  <- y[i-1]
        }
        results<-data.frame(cbind(y,lagy,x,lagx,z,e1,e2,nu))
        return(results) # drop first value
    }
    
    # create full panels
    dat<-do.call( rbind, replicate(units, gen.panel(1), simplify=FALSE ) )
    ## these functions might also work for x and u...
    #u<-as.vector(replicate(units,arima.sim(n=time,list(ar=phi))))
    #x<-as.vector(replicate(units,arima.sim(n=time,list(ar=rho))))+alpha*u
    
    # create plm object
    predictors<-data.table(dat,total)
    setkeyv(predictors,c('id','t'))
    data<-plm.data(predictors, c("id","t"))
    naive <- lm(y~x, data=data)
    naive_beta<- summary(naive)$coefficients[2,1]
    naive_se<- summary(naive)$coefficients[2,2]
    naive_t<- abs(naive_beta/naive_se)
    iv <- tsls(y  ~ x, instruments = ~ z, data=data)
    iv_beta<- summary(iv)$coefficients[2,1]
    iv_se<- summary(iv)$coefficients[2,2]
    iv_t<- abs(iv_beta/iv_se)
    lagid <- lm(y~lagx, data=data)
    lagid_beta<- summary(lagid)$coefficients[2,1]
    lagid_se<- summary(lagid)$coefficients[2,2]
    lagid_t<- abs(lagid_beta/lagid_se)
    #outputs <- c(phi,rho,alpha,beta,gamma,naive,iv,lagid)
    #outputs
    outputs <- cbind(phi,rho,alpha,beta,gamma,naive_beta,naive_se,naive_t,iv_beta,iv_se,iv_t,lagid_beta,lagid_se,lagid_t)
    return(outputs)
}


#lag.simulate<-Simulate(
#simulation(phi,rho,alpha,beta,gamma),
#expand.grid(
#phi<-c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
#rho<- .5, #define the list of values rho takes
#alpha<- .5, #define the list of values alpha takes.
#beta<- 2,
#gamma<- 1
#),
#nsim=1,
#trace=10
#)

phi_list<-c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
rho_list<-c(.5)
alpha_list<-c(-10,0,1,10)
alpha_list<-c(-10,0,1,10)
beta_list<-c(0,2)
gamma_list<-c(0,1,10)

nphi <- length(phi_list)
nrho <- length(rho_list)
nalpha <- length(alpha_list)
nbeta <- length(beta_list)
ngamma <- length(gamma_list)
numparams<- 15 #set the number of columns in the output (results as defined above) to which results are saved for each iteration
reps <- 100 #set the number of iteration
d <- data.frame(matrix(, nrow = nphi*nrho*nalpha*nbeta*ngamma*reps, ncol = numparams)) #d is the matrix that we save our results to.
i <- 0 #define i so that we can use it to save results to d[i,] for each iteration
output = NULL
for(phi in phi_list) {
    print(phi)
    for(rho in rho_list) {
        for(alpha in alpha_list) {
            print(alpha)
            for(beta in beta_list) {
                for(gamma in gamma_list) {
                    for(i3 in 1:reps) {
                        i <- i + 1
                        d[i,] <- cbind(i3,simulation(phi,rho,alpha,beta,gamma)) #for each iteration i3, which ranges from 0 through 1000, we save the results computed from model(n,rho_y,rho_x).
                    }
                }
            }
        }
    }
}

sim_results<-data.frame(d)
sim_results
names(sim_results) <- c("n","phi","rho","alpha","beta","gamma","naive_beta","naive_se","naive_t","iv_beta","iv_se","iv_t","lagid_beta","lagid_se","lagid_t")

write.table(sim_results, file ="sim_results/sim_results_iv.txt")

###########################################################################################################################################################################################
#Set gamma at 10
#plot
sim_results<-read.table("sim_results/sim_results_iv.txt", header=TRUE)
attach(sim_results)
sim_results_backup <- sim_results

sim_results <- sim_results_backup

sim_results1 <- sim_results[sim_results$alpha==0 & sim_results$beta==2 & sim_results$rho==.5 & sim_results$gamma==10, ]
sim_results2 <- sim_results[sim_results$alpha==1 & sim_results$beta==2 & sim_results$rho==.5 & sim_results$gamma==10, ]
sim_results3 <- sim_results[sim_results$alpha==10 & sim_results$beta==2 & sim_results$rho==.5 & sim_results$gamma==10, ]
sim_results4 <- sim_results[sim_results$alpha==-10 & sim_results$beta==2 & sim_results$rho==.5 & sim_results$gamma==10, ]

sim_results5 <- sim_results[sim_results$alpha==0 & sim_results$beta==0 & sim_results$rho==.5 & sim_results$gamma==10, ]
sim_results6 <- sim_results[sim_results$alpha==1 & sim_results$beta==0 & sim_results$rho==.5 & sim_results$gamma==10, ]
sim_results7 <- sim_results[sim_results$alpha==10 & sim_results$beta==0 & sim_results$rho==.5 & sim_results$gamma==10, ]
sim_results8 <- sim_results[sim_results$alpha==-10 & sim_results$beta==0 & sim_results$rho==.5 & sim_results$gamma==10, ]



attach(sim_results)
sim_results_backup <- sim_results
sim_results <- sim_results_backup

#prob of type2 error
sim_results1$iv_t_sig <- ifelse(sim_results1$iv_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results2$iv_t_sig <- ifelse(sim_results2$iv_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results3$iv_t_sig <- ifelse(sim_results3$iv_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results4$iv_t_sig <- ifelse(sim_results4$iv_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results1$naive_t_sig <- ifelse(sim_results1$naive_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results2$naive_t_sig <- ifelse(sim_results2$naive_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results3$naive_t_sig <- ifelse(sim_results3$naive_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results4$naive_t_sig <- ifelse(sim_results4$naive_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results1$lagid_t_sig <- ifelse(sim_results1$lagid_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results2$lagid_t_sig <- ifelse(sim_results2$lagid_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results3$lagid_t_sig <- ifelse(sim_results3$lagid_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results4$lagid_t_sig <- ifelse(sim_results4$lagid_t>=abs(qt(0.025,df=(5000-2))), 0, 1)

#prob of type 1 error
sim_results5$iv_t_sig <- ifelse(sim_results5$iv_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results6$iv_t_sig <- ifelse(sim_results6$iv_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results7$iv_t_sig <- ifelse(sim_results7$iv_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results8$iv_t_sig <- ifelse(sim_results8$iv_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results5$naive_t_sig <- ifelse(sim_results5$naive_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results6$naive_t_sig <- ifelse(sim_results6$naive_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results7$naive_t_sig <- ifelse(sim_results7$naive_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results8$naive_t_sig <- ifelse(sim_results8$naive_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results5$lagid_t_sig <- ifelse(sim_results5$lagid_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results6$lagid_t_sig <- ifelse(sim_results6$lagid_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results7$lagid_t_sig <- ifelse(sim_results7$lagid_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results8$lagid_t_sig <- ifelse(sim_results8$lagid_t>=abs(qt(0.025,df=(5000-2))), 1, 0)

ols.iv1 <- c(tapply(sim_results1$iv_beta, INDEX = sim_results1$phi, mean))#compute the mean of coefficients from the iv model
ols.naive1 <- c(tapply(sim_results1$naive_beta, INDEX = sim_results1$phi, mean))#compute the mean of naive estimates
ols.lagid1 <- c(tapply(sim_results1$lagid_beta, INDEX = sim_results1$phi, mean))#compute the mean of coefficients from the lag model

ols.iv2 <- c(tapply(sim_results2$iv_beta, INDEX = sim_results2$phi, mean))#compute the mean of coefficients from the iv model
ols.naive2 <- c(tapply(sim_results2$naive_beta, INDEX = sim_results2$phi, mean))#compute the mean of naive estimates
ols.lagid2 <- c(tapply(sim_results2$lagid_beta, INDEX = sim_results2$phi, mean))#compute the mean of coefficients from the lag model

ols.iv3 <- c(tapply(sim_results3$iv_beta, INDEX = sim_results3$phi, mean))#compute the mean of coefficients from the iv model
ols.naive3 <- c(tapply(sim_results3$naive_beta, INDEX = sim_results3$phi, mean))#compute the mean of naive estimates
ols.lagid3 <- c(tapply(sim_results3$lagid_beta, INDEX = sim_results3$phi, mean))#compute the mean of coefficients from the lag model

ols.iv4 <- c(tapply(sim_results4$iv_beta, INDEX = sim_results4$phi, mean))#compute the mean of coefficients from the iv model
ols.naive4 <- c(tapply(sim_results4$naive_beta, INDEX = sim_results4$phi, mean))#compute the mean of naive estimates
ols.lagid4 <- c(tapply(sim_results4$lagid_beta, INDEX = sim_results4$phi, mean))#compute the mean of coefficients from the lag model


ols.iv.sig1 <- c(tapply(sim_results1$iv_t_sig, INDEX = sim_results1$phi, mean))
ols.naive.sig1 <- c(tapply(sim_results1$naive_t_sig, INDEX = sim_results1$phi, mean))
ols.lagid.sig1 <- c(tapply(sim_results1$lagid_t_sig, INDEX = sim_results1$phi, mean))

ols.iv.sig2 <- c(tapply(sim_results2$iv_t_sig, INDEX = sim_results2$phi, mean))
ols.naive.sig2 <- c(tapply(sim_results2$naive_t_sig, INDEX = sim_results2$phi, mean))
ols.lagid.sig2 <- c(tapply(sim_results2$lagid_t_sig, INDEX = sim_results2$phi, mean))

ols.iv.sig3 <- c(tapply(sim_results3$iv_t_sig, INDEX = sim_results3$phi, mean))
ols.naive.sig3 <- c(tapply(sim_results3$naive_t_sig, INDEX = sim_results3$phi, mean))
ols.lagid.sig3 <- c(tapply(sim_results3$lagid_t_sig, INDEX = sim_results3$phi, mean))

ols.iv.sig4 <- c(tapply(sim_results4$iv_t_sig, INDEX = sim_results4$phi, mean))
ols.naive.sig4 <- c(tapply(sim_results4$naive_t_sig, INDEX = sim_results4$phi, mean))
ols.lagid.sig4 <- c(tapply(sim_results4$lagid_t_sig, INDEX = sim_results4$phi, mean))

ols.iv.sig5 <- c(tapply(sim_results5$iv_t_sig, INDEX = sim_results5$phi, mean))
ols.naive.sig5 <- c(tapply(sim_results5$naive_t_sig, INDEX = sim_results5$phi, mean))
ols.lagid.sig5 <- c(tapply(sim_results5$lagid_t_sig, INDEX = sim_results5$phi, mean))

ols.iv.sig6 <- c(tapply(sim_results6$iv_t_sig, INDEX = sim_results6$phi, mean))
ols.naive.sig6 <- c(tapply(sim_results6$naive_t_sig, INDEX = sim_results6$phi, mean))
ols.lagid.sig6  <- c(tapply(sim_results6$lagid_t_sig, INDEX = sim_results6$phi, mean))

ols.iv.sig7 <- c(tapply(sim_results7$iv_t_sig, INDEX = sim_results7$phi, mean))
ols.naive.sig7 <- c(tapply(sim_results7$naive_t_sig, INDEX = sim_results7$phi, mean))
ols.lagid.sig7  <- c(tapply(sim_results7$lagid_t_sig, INDEX = sim_results7$phi, mean))

ols.iv.sig8 <- c(tapply(sim_results8$iv_t_sig, INDEX = sim_results8$phi, mean))
ols.naive.sig8 <- c(tapply(sim_results8$naive_t_sig, INDEX = sim_results8$phi, mean))
ols.lagid.sig8  <- c(tapply(sim_results8$lagid_t_sig, INDEX = sim_results8$phi, mean))



ols.iv.bias1 <- ols.iv1 - 2 #2 is the true parameter
ols.naive.bias1 <- ols.naive1 - 2
ols.lagid.bias1 <- ols.lagid1 - 1 #1 is rho beta / (1-alpha beta) where beta = 2, rho = .5, alpha = 0
ols.iv.rmse1 <- (sim_results1$iv_beta - 2)^2 #2 is the true parameter
ols.naive.rmse1 <- (sim_results1$naive_beta - 2)^2
ols.lagid.rmse1 <- (sim_results1$lagid_beta - 1)^2 #1 is rho beta / (1-alpha beta) where beta = 2, rho = .5, alpha = 0
ols.iv.rmse1 <- sqrt((rowsum(ols.iv.rmse1,sim_results1$phi))/100)
ols.naive.rmse1 <- sqrt((rowsum(ols.naive.rmse1,sim_results1$phi))/100)
ols.lagid.rmse1 <- sqrt((rowsum(ols.lagid.rmse1,sim_results1$phi))/100)
bias1 <- c(ols.iv.bias1,ols.naive.bias1,ols.lagid.bias1)
rmse1 <- c(ols.iv.rmse1,ols.naive.rmse1,ols.lagid.rmse1)


ols.iv.bias2 <- ols.iv2 - 2 #2 is the true parameter
ols.naive.bias2 <- ols.naive2 - 2
ols.lagid.bias2 <- ols.lagid2 - 1 # -1 is rho beta / (1-alpha beta) where beta = 2, rho = .5, alpha = 1
ols.iv.rmse2 <- (sim_results2$iv_beta - 2)^2 #2 is the true parameter
ols.naive.rmse2 <- (sim_results2$naive_beta - 2)^2
ols.lagid.rmse2 <- (sim_results2$lagid_beta - 1)^2 # -1 is rho beta / (1-alpha beta) where beta = 2, rho = .5, alpha = 1
ols.iv.rmse2 <- sqrt((rowsum(ols.iv.rmse2,sim_results2$phi))/100)
ols.naive.rmse2 <- sqrt((rowsum(ols.naive.rmse2,sim_results2$phi))/100)
ols.lagid.rmse2 <- sqrt((rowsum(ols.lagid.rmse2,sim_results2$phi))/100)
bias2 <- c(ols.iv.bias2,ols.naive.bias2,ols.lagid.bias2)
rmse2 <- c(ols.iv.rmse2,ols.naive.rmse2,ols.lagid.rmse2)


ols.iv.bias3 <- ols.iv3 - 2
ols.naive.bias3 <- ols.naive3 - 2
ols.lagid.bias3 <- ols.lagid3 - 1 #is rho beta / (1-alpha beta) where beta = 2, rho = .5, alpha = 10
ols.iv.rmse3 <- (sim_results3$iv_beta - 2)^2
ols.naive.rmse3 <- (sim_results3$naive_beta - 2)^2
ols.lagid.rmse3 <- (sim_results3$lagid_beta -1)^2 #is rho beta / (1-alpha beta) where beta = 2, rho = .5, alpha = 10
ols.iv.rmse3 <- sqrt((rowsum(ols.iv.rmse3,sim_results3$phi))/100)
ols.naive.rmse3 <- sqrt((rowsum(ols.naive.rmse3,sim_results3$phi))/100)
ols.lagid.rmse3 <- sqrt((rowsum(ols.lagid.rmse3,sim_results3$phi))/100)
bias3 <- c(ols.iv.bias3,ols.naive.bias3,ols.lagid.bias3)
rmse3 <- c(ols.iv.rmse3,ols.naive.rmse3,ols.lagid.rmse3)


ols.iv.bias4 <- ols.iv4 - 2
ols.naive.bias4 <- ols.naive4 - 2
ols.lagid.bias4 <- ols.lagid4 - 0.047619048 #is rho beta / (1-alpha beta) where beta = 2, rho = .5, alpha = -10
ols.iv.rmse4 <- (sim_results4$iv_beta - 2)^2
ols.naive.rmse4 <- (sim_results4$naive_beta - 2)^2
ols.lagid.rmse4 <- (sim_results4$lagid_beta - 0.047619048)^2 #is rho beta / (1-alpha beta) where beta = 2, rho = .5, alpha = -10
ols.iv.rmse4 <- sqrt((rowsum(ols.iv.rmse4,sim_results4$phi))/100)
ols.naive.rmse4 <- sqrt((rowsum(ols.naive.rmse4,sim_results4$phi))/100)
ols.lagid.rmse4 <- sqrt((rowsum(ols.lagid.rmse4,sim_results4$phi))/100)
bias4 <- c(ols.iv.bias4,ols.naive.bias4,ols.lagid.bias4)
rmse4 <- c(ols.iv.rmse4,ols.naive.rmse4,ols.lagid.rmse4)


phi.vals<-(rep(c(seq(0,0.9,by=0.1)), times=3))
model<-(rep(c(1, 2, 3), each=10))

b1 <- bias1[model<4]
m1 <- model
m1 <-factor(m1, labels=c("IV", "NAIVE", "LAGID"))
data1 <- data.frame(b1,m1,phi.vals)
rmse1 <- rmse1[model<4]
data.rmse1 <- data.frame(rmse1,m1,phi.vals)

b2 <- bias2[model<4]
m2 <- model
m2 <-factor(m2, labels=c("IV", "NAIVE", "LAGID"))
data2 <- data.frame(b2,m2,phi.vals)
rmse2 <- rmse2[model<4]
data.rmse2 <- data.frame(rmse2,m2,phi.vals)

b3 <- bias3[model<4]
m3 <- model
m3 <-factor(m3, labels=c("IV", "NAIVE", "LAGID"))
data3 <- data.frame(b3,m3,phi.vals)
rmse3 <- rmse3[model<4]
data.rmse3 <- data.frame(rmse3,m3,phi.vals)

b4 <- bias4[model<4]
m4 <- model
m4 <-factor(m4, labels=c("IV", "NAIVE", "LAGID"))
data4 <- data.frame(b4,m4,phi.vals)
rmse4 <- rmse4[model<4]
data.rmse4 <- data.frame(rmse4,m4,phi.vals)

t1 <- c(ols.iv.sig1,ols.naive.sig1,ols.lagid.sig1)
m1 <- model
m1 <-factor(m1, labels=c("IV", "NAIVE", "LAGID"))
t1 <- t1[model<5]
data.t1 <- data.frame(t1,m1,phi.vals)

t2 <- c(ols.iv.sig2,ols.naive.sig2,ols.lagid.sig2)
m2 <- model
m2 <-factor(m2, labels=c("IV", "NAIVE", "LAGID"))
t2 <- t2[model<5]
data.t2 <- data.frame(t2,m2,phi.vals)

t3 <- c(ols.iv.sig3,ols.naive.sig3,ols.lagid.sig3)
m3 <- model
m3 <-factor(m3, labels=c("IV", "NAIVE", "LAGID"))
t3 <- t3[model<5]
data.t3 <- data.frame(t3,m3,phi.vals)

t4 <- c(ols.iv.sig4,ols.naive.sig4,ols.lagid.sig4)
m4 <- model
m4 <-factor(m4, labels=c("IV", "NAIVE", "LAGID"))
t4 <- t4[model<5]
data.t4 <- data.frame(t4,m4,phi.vals)

t5 <- c(ols.iv.sig5,ols.naive.sig5,ols.lagid.sig5)
m5 <- model
m5 <-factor(m5, labels=c("IV", "NAIVE", "LAGID"))
t5 <- t5[model<5]
data.t5 <- data.frame(t5,m5,phi.vals)

t6 <- c(ols.iv.sig6,ols.naive.sig6,ols.lagid.sig6)
m6 <- model
m6 <-factor(m6, labels=c("IV", "NAIVE", "LAGID"))
t6 <- t6[model<6]
data.t6 <- data.frame(t6,m6,phi.vals)

t7 <- c(ols.iv.sig7,ols.naive.sig7,ols.lagid.sig7)
m7 <- model
m7 <-factor(m7, labels=c("IV", "NAIVE", "LAGID"))
t7 <- t7[model<7]
data.t7 <- data.frame(t7,m7,phi.vals)

t8 <- c(ols.iv.sig8,ols.naive.sig8,ols.lagid.sig8)
m8 <- model
m8 <-factor(m8, labels=c("IV", "NAIVE", "LAGID"))
t8 <- t8[model<8]
data.t8 <- data.frame(t8,m8,phi.vals)


colours<-c(IV="green4", NAIVE="black", LAGID="red")
linetype<-c(IV="dotted", NAIVE="dashed", LAGID="solid")

#bias in estimated coefficients
fmt <- function(){
    f <- function(x) as.character(round(x,1))
    f
}

d<-ggplot(data1, aes(x=phi.vals, y=b1, group=m1))
plot1<-d+labs(title=expression(paste("Bias, ",alpha,"=0, ",beta,"=2")))+geom_line(aes(group=factor(m1), colour=factor(m1), linetype=factor(m1)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Bias", limits=c(-.2,.2), breaks=seq(from=-.2, to=.2, by=.1), labels=fmt())+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic() 

d<-ggplot(data2, aes(x=phi.vals, y=b2, group=m2))
plot2<-d+labs(title=expression(paste("Bias, ",alpha,"=1, ",beta,"=2")))+geom_line(aes(group=factor(m2), colour=factor(m2), linetype=factor(m2)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Bias", limits=c(-3,3), breaks=seq(from=-3, to=3, by=1))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data3, aes(x=phi.vals, y=b3, group=m3))
plot3<-d+labs(title=expression(paste("Bias, ",alpha,"=10, ",beta,"=2")))+geom_line(aes(group=factor(m3), colour=factor(m3), linetype=factor(m3)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Bias", limits=c(-3,3), breaks=seq(from=-3, to=3, by=1))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data4, aes(x=phi.vals, y=b4, group=m4))
plot4<-d+labs(title=expression(paste("Bias, ",alpha,"=-10, ",beta,"=2")))+geom_line(aes(group=factor(m4), colour=factor(m4), linetype=factor(m4)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Bias", limits=c(-3,3), breaks=seq(from=-3, to=3, by=1))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data.rmse1, aes(x=phi.vals, y=rmse1, group=m1))
plot5<-d+labs(title=expression(paste("RMSE, ",alpha,"=0, ",beta,"=2")))+geom_line(aes(group=factor(m1), colour=factor(m1), linetype=factor(m1)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("RMSE", limits=c(0,.3), breaks=seq(from=0, to=.3, by=0.1))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data.rmse2, aes(x=phi.vals, y=rmse2, group=m2))
plot6<-d+labs(title=expression(paste("RMSE, ",alpha,"=1, ",beta,"=2")))+geom_line(aes(group=factor(m2), colour=factor(m2), linetype=factor(m2)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("RMSE", limits=c(0,3), breaks=seq(from=0, to=3, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data.rmse3, aes(x=phi.vals, y=rmse3, group=m3))
plot7<-d+labs(title=expression(paste("RMSE, ",alpha,"=10, ",beta,"=2")))+geom_line(aes(group=factor(m3), colour=factor(m3), linetype=factor(m3)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("RMSE", limits=c(0,3), breaks=seq(from=0, to=3, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data.rmse4, aes(x=phi.vals, y=rmse4, group=m4))
plot8<-d+labs(title=expression(paste("RMSE, ",alpha,"=-10, ",beta,"=2")))+geom_line(aes(group=factor(m4), colour=factor(m4), linetype=factor(m4)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("RMSE", limits=c(0,3), breaks=seq(from=0, to=3, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t5, aes(x=phi.vals, y=t5, group=m5))
plot9<-f+labs(title=expression(paste("Type-1 Error, ",alpha,"=0, ",beta,"=0")))+geom_line(aes(group=factor(m5), colour=factor(m5), linetype=factor(m5)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t6, aes(x=phi.vals, y=t6, group=m6))
plot10<-f+labs(title=expression(paste("Type-1 Error, ",alpha,"=1, ",beta,"=0")))+geom_line(aes(group=factor(m6), colour=factor(m6), linetype=factor(m6)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t7, aes(x=phi.vals, y=t7, group=m7))
plot11<-f+labs(title=expression(paste("Type-1 Error, ",alpha,"=10, ",beta,"=0")))+geom_line(aes(group=factor(m7), colour=factor(m7), linetype=factor(m7)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t8, aes(x=phi.vals, y=t8, group=m8))
plot12<-f+labs(title=expression(paste("Type-1 Error, ",alpha,"=-10, ",beta,"=0")))+geom_line(aes(group=factor(m8), colour=factor(m8), linetype=factor(m8)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t1, aes(x=phi.vals, y=t1, group=m8))
plot13<-f+labs(title=expression(paste("Type-2 Error, ",alpha,"=0, ",beta,"=2")))+geom_line(aes(group=factor(m1), colour=factor(m1), linetype=factor(m1)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t2, aes(x=phi.vals, y=t2, group=m8))
plot14<-f+labs(title=expression(paste("Type-2 Error, ",alpha,"=1, ",beta,"=2")))+geom_line(aes(group=factor(m2), colour=factor(m2), linetype=factor(m2)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t3, aes(x=phi.vals, y=t3, group=m8))
plot15<-f+labs(title=expression(paste("Type-2 Error, ",alpha,"=10, ",beta,"=2")))+geom_line(aes(group=factor(m3), colour=factor(m3), linetype=factor(m3)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t4, aes(x=phi.vals, y=t4, group=m8))
plot16<-f+labs(title=expression(paste("Type-2 Error, ",alpha,"=-10, ",beta,"=2")))+geom_line(aes(group=factor(m4), colour=factor(m4), linetype=factor(m4)), size=.8) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

#### combine plots, one legend
library(ggplot2)
library(gridExtra)

grid_arrange_shared_legend <- function(...) {
  plots <- list(...)
  g <- ggplotGrob(plots[[1]] + theme(legend.position="bottom"))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  grid.arrange(arrangeGrob(grobs=lapply(plots, function(x)
      x + theme(legend.position="none")),ncol = 3),
    legend,
    ncol = 1,
    heights = unit.c(unit(1, "npc") - lheight, lheight))
}

grid_arrange_shared_legend(plot2, plot6, plot10, plot3, plot7, plot11)
plotted<-recordPlot()
pdf("figures/fig9.pdf", width=9, height=8)
replayPlot(plotted)
dev.off()

###########################################################################################################################################################################################
#Set gamma at 0
#plot
sim_results<-read.table("sim_results/sim_results_iv.txt", header=TRUE)
attach(sim_results)
sim_results_backup <- sim_results

sim_results <- sim_results_backup

sim_results1 <- sim_results[sim_results$alpha==0 & sim_results$beta==2 & sim_results$rho==.5 & sim_results$gamma==0, ]
sim_results2 <- sim_results[sim_results$alpha==1 & sim_results$beta==2 & sim_results$rho==.5 & sim_results$gamma==0, ]
sim_results3 <- sim_results[sim_results$alpha==10 & sim_results$beta==2 & sim_results$rho==.5 & sim_results$gamma==0, ]
sim_results4 <- sim_results[sim_results$alpha==-10 & sim_results$beta==2 & sim_results$rho==.5 & sim_results$gamma==0, ]

sim_results5 <- sim_results[sim_results$alpha==0 & sim_results$beta==0 & sim_results$rho==.5 & sim_results$gamma==0, ]
sim_results6 <- sim_results[sim_results$alpha==1 & sim_results$beta==0 & sim_results$rho==.5 & sim_results$gamma==0, ]
sim_results7 <- sim_results[sim_results$alpha==10 & sim_results$beta==0 & sim_results$rho==.5 & sim_results$gamma==0, ]
sim_results8 <- sim_results[sim_results$alpha==-10 & sim_results$beta==0 & sim_results$rho==.5 & sim_results$gamma==0, ]



attach(sim_results)
sim_results_backup <- sim_results
sim_results <- sim_results_backup

#prob of type2 error
sim_results1$naive_t_sig <- ifelse(sim_results1$naive_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results2$naive_t_sig <- ifelse(sim_results2$naive_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results3$naive_t_sig <- ifelse(sim_results3$naive_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results4$naive_t_sig <- ifelse(sim_results4$naive_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results1$lagid_t_sig <- ifelse(sim_results1$lagid_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results2$lagid_t_sig <- ifelse(sim_results2$lagid_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results3$lagid_t_sig <- ifelse(sim_results3$lagid_t>=abs(qt(0.025,df=(5000-2))), 0, 1)
sim_results4$lagid_t_sig <- ifelse(sim_results4$lagid_t>=abs(qt(0.025,df=(5000-2))), 0, 1)

#prob of type 1 error
sim_results5$naive_t_sig <- ifelse(sim_results5$naive_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results6$naive_t_sig <- ifelse(sim_results6$naive_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results7$naive_t_sig <- ifelse(sim_results7$naive_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results8$naive_t_sig <- ifelse(sim_results8$naive_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results5$lagid_t_sig <- ifelse(sim_results5$lagid_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results6$lagid_t_sig <- ifelse(sim_results6$lagid_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results7$lagid_t_sig <- ifelse(sim_results7$lagid_t>=abs(qt(0.025,df=(5000-2))), 1, 0)
sim_results8$lagid_t_sig <- ifelse(sim_results8$lagid_t>=abs(qt(0.025,df=(5000-2))), 1, 0)

ols.naive1 <- c(tapply(sim_results1$naive_beta, INDEX = sim_results1$phi, mean))#compute the mean of naive estimates
ols.lagid1 <- c(tapply(sim_results1$lagid_beta, INDEX = sim_results1$phi, mean))#compute the mean of coefficients from the lag model

ols.naive2 <- c(tapply(sim_results2$naive_beta, INDEX = sim_results2$phi, mean))#compute the mean of naive estimates
ols.lagid2 <- c(tapply(sim_results2$lagid_beta, INDEX = sim_results2$phi, mean))#compute the mean of coefficients from the lag model

ols.naive3 <- c(tapply(sim_results3$naive_beta, INDEX = sim_results3$phi, mean))#compute the mean of naive estimates
ols.lagid3 <- c(tapply(sim_results3$lagid_beta, INDEX = sim_results3$phi, mean))#compute the mean of coefficients from the lag model

ols.naive4 <- c(tapply(sim_results4$naive_beta, INDEX = sim_results4$phi, mean))#compute the mean of naive estimates
ols.lagid4 <- c(tapply(sim_results4$lagid_beta, INDEX = sim_results4$phi, mean))#compute the mean of coefficients from the lag model


ols.naive.sig1 <- c(tapply(sim_results1$naive_t_sig, INDEX = sim_results1$phi, mean))
ols.lagid.sig1 <- c(tapply(sim_results1$lagid_t_sig, INDEX = sim_results1$phi, mean))

ols.naive.sig2 <- c(tapply(sim_results2$naive_t_sig, INDEX = sim_results2$phi, mean))
ols.lagid.sig2 <- c(tapply(sim_results2$lagid_t_sig, INDEX = sim_results2$phi, mean))

ols.naive.sig3 <- c(tapply(sim_results3$naive_t_sig, INDEX = sim_results3$phi, mean))
ols.lagid.sig3 <- c(tapply(sim_results3$lagid_t_sig, INDEX = sim_results3$phi, mean))

ols.naive.sig4 <- c(tapply(sim_results4$naive_t_sig, INDEX = sim_results4$phi, mean))
ols.lagid.sig4 <- c(tapply(sim_results4$lagid_t_sig, INDEX = sim_results4$phi, mean))

ols.naive.sig5 <- c(tapply(sim_results5$naive_t_sig, INDEX = sim_results5$phi, mean))
ols.lagid.sig5 <- c(tapply(sim_results5$lagid_t_sig, INDEX = sim_results5$phi, mean))

ols.naive.sig6 <- c(tapply(sim_results6$naive_t_sig, INDEX = sim_results6$phi, mean))
ols.lagid.sig6  <- c(tapply(sim_results6$lagid_t_sig, INDEX = sim_results6$phi, mean))

ols.naive.sig7 <- c(tapply(sim_results7$naive_t_sig, INDEX = sim_results7$phi, mean))
ols.lagid.sig7  <- c(tapply(sim_results7$lagid_t_sig, INDEX = sim_results7$phi, mean))

ols.naive.sig8 <- c(tapply(sim_results8$naive_t_sig, INDEX = sim_results8$phi, mean))
ols.lagid.sig8  <- c(tapply(sim_results8$lagid_t_sig, INDEX = sim_results8$phi, mean))



ols.naive.bias1 <- ols.naive1 - 2
ols.lagid.bias1 <- ols.lagid1 - 1 #1 is beta*rho where beta is set at 2 and rho at .5.
ols.naive.rmse1 <- (sim_results1$naive_beta - 2)^2
ols.lagid.rmse1 <- (sim_results1$lagid_beta - 1)^2 #1 is beta*rho where beta is set at 2 and rho at .5.
ols.naive.rmse1 <- sqrt((rowsum(ols.naive.rmse1,sim_results1$phi))/100)
ols.lagid.rmse1 <- sqrt((rowsum(ols.lagid.rmse1,sim_results1$phi))/100)
bias1 <- c(ols.naive.bias1,ols.lagid.bias1)
rmse1 <- c(ols.naive.rmse1,ols.lagid.rmse1)


ols.naive.bias2 <- ols.naive2 - 2
ols.lagid.bias2 <- ols.lagid2 - 1 #1 is beta*rho where beta is set at 2 and rho at 0.
ols.naive.rmse2 <- (sim_results2$naive_beta - 2)^2
ols.lagid.rmse2 <- (sim_results2$lagid_beta - 1)^2 #1 is beta*rho where beta is set at 2 and rho at 0.
ols.naive.rmse2 <- sqrt((rowsum(ols.naive.rmse2,sim_results2$phi))/100)
ols.lagid.rmse2 <- sqrt((rowsum(ols.lagid.rmse2,sim_results2$phi))/100)
bias2 <- c(ols.naive.bias2,ols.lagid.bias2)
rmse2 <- c(ols.naive.rmse2,ols.lagid.rmse2)


ols.naive.bias3 <- ols.naive3 - 2
ols.lagid.bias3 <- ols.lagid3 - 1 #1 is beta*rho where beta is set at 2 and rho at 0.5.
ols.naive.rmse3 <- (sim_results3$naive_beta - 2)^2
ols.lagid.rmse3 <- (sim_results3$lagid_beta - 1)^2
ols.naive.rmse3 <- sqrt((rowsum(ols.naive.rmse3,sim_results3$phi))/100)
ols.lagid.rmse3 <- sqrt((rowsum(ols.lagid.rmse3,sim_results3$phi))/100)
bias3 <- c(ols.naive.bias3,ols.lagid.bias3)
rmse3 <- c(ols.naive.rmse3,ols.lagid.rmse3)


ols.naive.bias4 <- ols.naive4 - 2
ols.lagid.bias4 <- ols.lagid4 - 1 #1 is beta*rho where beta is set at 2 and rho at 0.5.
ols.naive.rmse4 <- (sim_results4$naive_beta - 2)^2
ols.lagid.rmse4 <- (sim_results4$lagid_beta - 1)^2
ols.naive.rmse4 <- sqrt((rowsum(ols.naive.rmse4,sim_results4$phi))/100)
ols.lagid.rmse4 <- sqrt((rowsum(ols.lagid.rmse4,sim_results4$phi))/100)
bias4 <- c(ols.naive.bias4,ols.lagid.bias4)
rmse4 <- c(ols.naive.rmse4,ols.lagid.rmse4)


phi.vals<-(rep(c(seq(0,0.9,by=0.1)), times=2))
model<-(rep(c(1, 2), each=10))

b1 <- bias1[model<3]
m1 <- model
m1 <-factor(m1, labels=c("NAIVE", "LAGID"))
data1 <- data.frame(b1,m1,phi.vals)
rmse1 <- rmse1[model<3]
data.rmse1 <- data.frame(rmse1,m1,phi.vals)

b2 <- bias2[model<3]
m2 <- model
m2 <-factor(m2, labels=c("NAIVE", "LAGID"))
data2 <- data.frame(b2,m2,phi.vals)
rmse2 <- rmse2[model<3]
data.rmse2 <- data.frame(rmse2,m2,phi.vals)

b3 <- bias3[model<3]
m3 <- model
m3 <-factor(m3, labels=c("NAIVE", "LAGID"))
data3 <- data.frame(b3,m3,phi.vals)
rmse3 <- rmse3[model<3]
data.rmse3 <- data.frame(rmse3,m3,phi.vals)

b4 <- bias4[model<3]
m4 <- model
m4 <-factor(m4, labels=c("NAIVE", "LAGID"))
data4 <- data.frame(b4,m4,phi.vals)
rmse4 <- rmse4[model<3]
data.rmse4 <- data.frame(rmse4,m4,phi.vals)

t1 <- c(ols.naive.sig1,ols.lagid.sig1)
m1 <- model
m1 <-factor(m1, labels=c("NAIVE", "LAGID"))
t1 <- t1[model<3]
data.t1 <- data.frame(t1,m1,phi.vals)

t2 <- c(ols.naive.sig2,ols.lagid.sig2)
m2 <- model
m2 <-factor(m2, labels=c("NAIVE", "LAGID"))
t2 <- t2[model<3]
data.t2 <- data.frame(t2,m2,phi.vals)

t3 <- c(ols.naive.sig3,ols.lagid.sig3)
m3 <- model
m3 <-factor(m3, labels=c("NAIVE", "LAGID"))
t3 <- t3[model<3]
data.t3 <- data.frame(t3,m3,phi.vals)

t4 <- c(ols.naive.sig4,ols.lagid.sig4)
m4 <- model
m4 <-factor(m4, labels=c("NAIVE", "LAGID"))
t4 <- t4[model<3]
data.t4 <- data.frame(t4,m4,phi.vals)

t5 <- c(ols.naive.sig5,ols.lagid.sig5)
m5 <- model
m5 <-factor(m5, labels=c("NAIVE", "LAGID"))
t5 <- t5[model<3]
data.t5 <- data.frame(t5,m5,phi.vals)

t6 <- c(ols.naive.sig6,ols.lagid.sig6)
m6 <- model
m6 <-factor(m6, labels=c("NAIVE", "LAGID"))
t6 <- t6[model<3]
data.t6 <- data.frame(t6,m6,phi.vals)

t7 <- c(ols.naive.sig7,ols.lagid.sig7)
m7 <- model
m7 <-factor(m7, labels=c("NAIVE", "LAGID"))
t7 <- t7[model<3]
data.t7 <- data.frame(t7,m7,phi.vals)

t8 <- c(ols.naive.sig8,ols.lagid.sig8)
m8 <- model
m8 <-factor(m8, labels=c("NAIVE", "LAGID"))
t8 <- t8[model<3]
data.t8 <- data.frame(t8,m8,phi.vals)


colours<-c(NAIVE="black", LAGID="red")
linetype<-c(NAIVE="dashed", LAGID="solid")


#bias in estimated coefficients
fmt <- function(){
  f <- function(x) as.character(round(x,1))
  f
}

d<-ggplot(data1, aes(x=phi.vals, y=b1, group=m1))
plot1<-d+labs(title=expression(paste("Bias, ",alpha,"=0, ",beta,"=2")))+geom_line(aes(group=factor(m1), colour=factor(m1), linetype=factor(m1)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Bias", limits=c(-3,3), breaks=seq(from=-3, to=3, by=1))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic() 

d<-ggplot(data2, aes(x=phi.vals, y=b2, group=m2))
plot2<-d+labs(title=expression(paste("Bias, ",alpha,"=1, ",beta,"=2")))+geom_line(aes(group=factor(m2), colour=factor(m2), linetype=factor(m2)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Bias", limits=c(-3,3), breaks=seq(from=-3, to=3, by=1))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data3, aes(x=phi.vals, y=b3, group=m3))
plot3<-d+labs(title=expression(paste("Bias, ",alpha,"=10, ",beta,"=2")))+geom_line(aes(group=factor(m3), colour=factor(m3), linetype=factor(m3)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Bias", limits=c(-3,3), breaks=seq(from=-3, to=3, by=1))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data4, aes(x=phi.vals, y=b4, group=m4))
plot4<-d+labs(title=expression(paste("Bias, ",alpha,"=-10, ",beta,"=2")))+geom_line(aes(group=factor(m4), colour=factor(m4), linetype=factor(m4)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Bias", limits=c(-3,3), breaks=seq(from=-3, to=3, by=1))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data.rmse1, aes(x=phi.vals, y=rmse1, group=m1))
plot5<-d+labs(title=expression(paste("RMSE, ",alpha,"=0, ",beta,"=2")))+geom_line(aes(group=factor(m1), colour=factor(m1), linetype=factor(m1)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("RMSE", limits=c(0,3), breaks=seq(from=0, to=3, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data.rmse2, aes(x=phi.vals, y=rmse2, group=m2))
plot6<-d+labs(title=expression(paste("RMSE, ",alpha,"=1, ",beta,"=2")))+geom_line(aes(group=factor(m2), colour=factor(m2), linetype=factor(m2)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("RMSE", limits=c(0,3), breaks=seq(from=0, to=3, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data.rmse3, aes(x=phi.vals, y=rmse3, group=m3))
plot7<-d+labs(title=expression(paste("RMSE, ",alpha,"=10, ",beta,"=2")))+geom_line(aes(group=factor(m3), colour=factor(m3), linetype=factor(m3)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("RMSE", limits=c(0,3), breaks=seq(from=0, to=3, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

d<-ggplot(data.rmse4, aes(x=phi.vals, y=rmse4, group=m4))
plot8<-d+labs(title=expression(paste("RMSE, ",alpha,"=-10, ",beta,"=2")))+geom_line(aes(group=factor(m4), colour=factor(m4), linetype=factor(m4)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("RMSE", limits=c(0,3), breaks=seq(from=0, to=3, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t5, aes(x=phi.vals, y=t5, group=m5))
plot9<-f+labs(title=expression(paste("Type-1 Error, ",alpha,"=0, ",beta,"=0")))+geom_line(aes(group=factor(m5), colour=factor(m5), linetype=factor(m5)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t6, aes(x=phi.vals, y=t6, group=m6))
plot10<-f+labs(title=expression(paste("Type-1 Error, ",alpha,"=1, ",beta,"=0")))+geom_line(aes(group=factor(m6), colour=factor(m6), linetype=factor(m6)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t7, aes(x=phi.vals, y=t7, group=m7))
plot11<-f+labs(title=expression(paste("Type-1 Error, ",alpha,"=10, ",beta,"=0")))+geom_line(aes(group=factor(m7), colour=factor(m7), linetype=factor(m7)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t8, aes(x=phi.vals, y=t8, group=m8))
plot12<-f+labs(title=expression(paste("Type-1 Error, ",alpha,"=-10, ",beta,"=0")))+geom_line(aes(group=factor(m8), colour=factor(m8), linetype=factor(m8)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t1, aes(x=phi.vals, y=t1, group=m8))
plot13<-f+labs(title=expression(paste("Type-2 Error, ",alpha,"=0, ",beta,"=2")))+geom_line(aes(group=factor(m1), colour=factor(m1), linetype=factor(m1)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t2, aes(x=phi.vals, y=t2, group=m8))
plot14<-f+labs(title=expression(paste("Type-2 Error, ",alpha,"=1, ",beta,"=2")))+geom_line(aes(group=factor(m2), colour=factor(m2), linetype=factor(m2)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t3, aes(x=phi.vals, y=t3, group=m8))
plot15<-f+labs(title=expression(paste("Type-2 Error, ",alpha,"=10, ",beta,"=2")))+geom_line(aes(group=factor(m3), colour=factor(m3), linetype=factor(m3)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

f<-ggplot(data.t4, aes(x=phi.vals, y=t4, group=m8))
plot16<-f+labs(title=expression(paste("Type-2 Error, ",alpha,"=-10, ",beta,"=2")))+geom_line(aes(group=factor(m4), colour=factor(m4), linetype=factor(m4)), size=.5) +scale_x_continuous(expression(paste(phi)), limits=c(0,0.9), labels=fmt())+scale_y_continuous("Frequency", limits=c(0,1), breaks=seq(from=0, to=1, by=0.5))+ scale_colour_manual(name="Method", values=colours)+scale_linetype_manual(name="Method", values=linetype)+geom_abline(intercept=0, slope=0, colour="darkgray") + theme_classic()

#### combine plots, one legend
library(ggplot2)
library(gridExtra)

grid_arrange_shared_legend <- function(...) {
  plots <- list(...)
  g <- ggplotGrob(plots[[1]] + theme(legend.position="bottom"))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  grid.arrange(
    do.call(arrangeGrob, lapply(plots, function(x)
      x + theme(legend.position="none"))),
    legend,
    ncol = 1,
    heights = unit.c(unit(1, "npc") - lheight, lheight))
}

grid_arrange_shared_legend(plot1, plot2, plot3, plot5, plot6, plot7, plot9, plot10, plot11)
plotted<-recordPlot()
pdf("figures/figA-4.pdf", width=8, height=9)
replayPlot(plotted)
dev.off()




